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Abstract 

Results are presented from numerical experiments aiming at the computation of 
stochastic phase-field models for phase transformations by coarse-graining molecular 
dynamics. The studied phase transformations occur between a solid crystal and a 
liquid. Nuclcation and growth, sometimes dendritic, of crystal grains in a sub-cooled 
liquid is determined by diffusion and convection of heat, on the macroscopic level, 
and by interface effects, where the width of the solid-liquid interface is on an atomic 
length-scale. Phase-field methods are widely used in the study of the continuum level 
time evolution of the phase transformations; they introduce an order parameter to 
distinguish between the phases. The dynamics of the order parameter is modelled by 
an Allen-Cahn equation and coupled to an energy equation, where the latent heat at 
the phase transition enters as a source term. Stochastic fluctuations are sometimes 
added in the coupled system of partial differential equations to introdxicc nuclcation 
and to get qualitatively correct behaviour of dendritic side-branching. In this report 
the possibility of computing some of the Allen-Cahn model functions from a microscale 
model is investigated. The microscopic model description of the material by stochas- 
tic, Smoluchowski, dynamics is considered given. A local average of contributions to 
the potential energy in the micro model is used to determine the local phase, and a 
stochastic phase-field model is computed by coarse-graining the molecular dynamics. 
Molecular dynamics simulations on a two phase system at the melting point are used 
to compute a double-well reaction term in the Allen-Cahn equation and a diffusion 
matrix describing the noise in the coarse-grained phase-field. 

This work was supported by the Swedish Foundation for Strategic Research grant A3 02:123, 
"Mathematical theory and simulation tools for phase transformations is materials". 



1 Introduction 



Phase-field methods are widely used for modelling phase transformations in materials on 
the continuum level and exist in many different versions for different applications. In this 



report the considered phase transformation occurs in a single component system with a 
solid and a liquid phase. 

The phase-field model of solidification studied here is a coupled system of partial differential 
equations for the temperature, T, and a phase-field, 0, which is an order parameter used to 
distinguish between the solid and the liquid subdomains. Two different values, (j)s and 4>i, are 
equilibrium values of the phase-field in solid and liquid respectively. The phase-field varies 
continuously between the two values and the interface between solid and liquid, at a time t, 
is defined as a level surface of the phase-field; for example {x € M'' : (j){x, t) = 0.5{(j>s + '/'/)}• 
From a computational point of view the implicit definition of the phases in the phase-field 
method, as in the level set method [5J [T^], is an advantage over sharp interface methods, 
since it avoids the explicit tracking of the interface. A local change of the phase-field from 
to (j)s in a subdomain translates into solidification of that region with a corresponding 
release of latent heat and the reverse change from (ps to (pi means melting which requires 
energy. The release or absorption of latent heat is modelled as a continuous function of (j) 
so that the energy released when a unit volume solidifies is L{g{(j)i) — g{4>s)), where L is the 
latent heat and g{(j)) is a model function, monotone with g{4>s) = 0, g{4>i) = 1, 9'{4>s) = 0, 
and g'{4>i) = 0. Then the energy equation for a unit volume becomes a heat equation with 
a source term 

^(cvr + Lg(0)) = V-(AVT), 

where cy is the heat capacity at constant volume and A is the thermal conductivity. Here, 
and in the following, the usual notation for differentiation with respect to the spatial vari- 
ables is applied, with V and V- denoting the gradient and the divergence respectively. The 
phase-field, and the related model function g, are exceptional in the energy equation in the 
sense that, while all the other quantities are standard physical quantities on the macroscopic 
level, the phase-field need not be associated with a measurable quantity. A phenomenolog- 
ical model of the phase change is given by the energy equation coupled to the AUen-Cahn 
equation 

^ = V • (fciV0) - k2 [f\4>) + .g'(0)A:3(TM - T)) (1) 

for the time evolution of the phase-field; here Tm denotes the melting point, fci, fc2, and 
^3, are positive model parameters (fci may be an anisotropic matrix introducing directional 
dependence on the growth of the solid) , and the model function / is a double well potential 
with minima at <f>s and (j)i . Standard examples of the model functions are 

when 0<; = — 1 and (pi = \. By construction of the model functions, the reaction term in the 
AUen-Cahn equation vanishes where (f> — (f>s or (j) — independently of the temperature. 
Since the diffusion term is zero for any constant function the two constant phase-fields 
(j> = (j>s and (j) = (j>i are stationary solutions to the Allen-Cahn equation for all temperatures. 
This means, for example, that nucleation of solid in a region of subcooled liquid can not 
occur in a phase-field modelled by the deterministic Allen-Cahn equation above. The effect 
of nucleation can be introduced in the model by adding a noise term in the Allen-Cahn 
equation, giving a stochastic partial differential equation. Simulation of dendrite growth in 
an subcooled liquid is another example where the deterministic system is inadequate; its 
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solutions fail to develop the side branches seen to form in real dendrites as the tips grow. 
Stochastic phase-field models where noise is added to either one, or both, of the AUen-Cahn 
equation and the energy equation are used to include the effect of side branching; see for 
example [2]- 

The present report contains the results from numerical experiments on a method presented 
and analysed in |T3j and the rest of this introduction summarises the ideas from jl4j needed 
here. That report takes the stochastic phase-field model 

^(cvr+L5(^)) = V-(AVr), (2a) 

^ = V • (fci V0) - fc2 (/'('^) + g\cf>)h{TM ~ T)) + noise, (2b) 

as its starting point and asks whether it is possible to obtain the model functions and 
parameters, /(0), g{(l>)j ^2) ^3, and the noise, from computations on a microscale model. 
To answer this question the phase-field, 0, must be defined in terms of quantities computable 
on the microscale. The microscopic model used for this purpose is a molecular dynamics 
model of N particles in a microscopic domain D in M'^ where the motion of the particles is 
given by the Smoluchowski dynamics; see for example |5]. Thus, with X* e denoting 
the positions of all particles in the system at the time t and g M'^ the position of particle 
i, the dynamics are given by the Ito stochastic differential equations 

dXl = -VxViX*) dt + y/2kBT dWl i= l,2,...,iV, (3) 

where U is the total potential energy of the system, V x. denotes the gradient with respect 
to the position of particle i, k-g, is the Boltzmann constant, and Wi = (W^.i, Wi^2, W^i.s)"'" 
are independent three dimensional Brownian motions, with independent components. The 
macroscopic temperature, T, is a constant input parameter in the microscopic model. We 
may identify the latent heat, in the macroscopic model, with the difference in total potential 
energy per unit volume of the liquid and the solid at the melting point, in the microscopic 
model. The idea is then to let the local contributions to the total potential energy define the 
phase variable. Since the potential energy decreases with the temperature even in a single 
phase system the equilibrium values of such a phase-field, m, unlike those of 0, depend on 
the temperature; see Figure [T] Assuming that in pure solid or pure liquid the phase-field. 



Subcooled Liquid Liquid 



Latent lieat, L 



Solid Superlieated Solid 



-*■ T 



Figure 1: Schematic picture of m(T) for a pure liquid (top curve) and a pure solid (bottom 
curve) and the latent heat as the jump in m at a phase transition. 

m, varies slowly, compared to the latent heat release, with the temperature close to the 
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melting point, the energy equation becomes 

d 



{cvT + m)^V -{XVT), 



where cy and A are approximately the same as in ( 2a I for T !v Tj 



M- 



For a model where the total potential energy of the system can be naturally split into a sum 
of contributions arising from the interaction of individual atoms with their environment, 



N 



U{X) = Y.m,{X), (4) 

i=l 

phase-fields can be introduced on the micro level by localised averages of these contributions; 
a given configuration X defines a phase-field m( • ; X) : D ^ M. through 

N 

m(a;;X) = ^m,(X),7(x-X,), (5) 

where the choice of moUifier, rj, determines the spatial smoothness of the phase-field. If, for 
example, the potential energy is defined entirely by pairwise interactions 

N N 

as is common in simple molecular dynamics models, it is natural to let 

N 



be particle z's contribution to the total potential energy. 

With the definition ^ of the potential energy phase-field, m, and with the microscopic 
system defined by ([3| and Q , Ito's formula gives a stochastic differential equation 

N 3 

dm{x;X*) = a{x;X*) dt + J^^l^jA^'^^*) dWl^, (6) 

j=i k=i 

for m evaluated in a point x ^ D. The drift, a(x] •), and the diffusions, [3j^k{x\ •), are explic- 
itly known functions expressed in terms of the m^is, the moUifier, ry, and their derivatives 
up to second order. While m by definition is a continuous field it is still an atomic scale 
quantity since it is defined in terms the particle positions X*. A macroscopic phase-field, 
similar to in (|2]), must lose both the dependence on the particle positions, X*, and the 
explicit dependence on the microscale space variable x. To achieve this, a coarse-grained 
approximation mcg{x) of m{x) is introduced as a solution of a stochastic differential equa- 
tion 

M 

dml^ix) = a(m*g)(x) dt + ^6,(m*g)(x) dW^ , (7) 
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where the independent Wiener processes Wj, j — 1,2, . . . , Af ^ N, also are independent 
of the Wiener processes Wi in the micro model. Here the drift and diffusion coefficient 
functions, a(m*g) and bj^m^.^), may depend on more information about the coarse-grained 
phase-field than just the point value; compare the stochastic AUen-Cahn equation (2b), 
where the diffusion term in the drift contains second derivatives of the phase-field. 



The choice of the coarse-grained drift and diffusion functions proceeds in two steps: first, 
finding a general form the coarse-grained equation where the drift and diffusion coefficient 
functions, defined as time averaged expected values of the microscopic drift and diffusions 
over simulation paths, still depend on the micro scale space variable, x; second, expressing 
the X dependent coarse-grained drift and diffusion coefficients by drift and diffusion func- 
tions depending only on the phase- field meg, using that Wcg is a smooth monotone function 
of X in the interface. 

In the first step, a coarse-grained stochastic differential equation 

M 



dm*g(x) = a{x) dt + Y^bj{x) dW* , 



is introduced by defining the drift 



a{x) = -E 



a{x;X*)dt 







X e D, 



(8a) 



and choosing a diffusion matrix that fulfil 

M r „r 3 

Yb,{x)-b,{x') = -E / 

7 = 1 r .7 = 1 k = l 



*)f3,,k{x';X*)dt 



X^=Xr. 



X, x' e D, (8b) 



for some fixed, deterministic, initial conditions X^ = Xq. The initial condition for the 
coarse-grained phase- field is m^g = to(-;^o)- This particular coarse-graining is motivated 
by the argument that the coarse-grained model will be used to compute properties on the 
form E [y(m(-; X^))] , where y : D ^ M is a smooth function and T > is a fixed final 
time. The optimal coarse-grained model is the one that minimises the error in the expected 
value; using the conditional expected values u(/i,t) = E[y(m^) | m*g = /i], this error can 
be expressed as 

E[yim{-,X^))]-E[yiml)] 



LHD) 



iu'{m{-,X'),t) , a{--X')-a{-))^ 
- / /u"(m(.;X*),i) , ^^(/3,-fe® )(.,.; X*)- ^(5, »&,)(.,.; 

"'O \ j = lfc=l i = \ 



dt 



L^{DxD) 



where (E) denotes the tensor product (bj (Ebj){x, x') = bj{x)bj{x'), and u' and u" denote the 
first and second variations of u(/i, t) with respect to fi. Assuming that u' can be expanded 
in powers of a — a, the choice (8a I cancels the leading term in the error associated with u' . 



Similarly, (8b I corresponds to cancelling the dominating term in the expansion of u". 
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In a practical computation the functions a and f3j can only be evaluated in a discrete set 



of points Dk = {x^, ■ ■ ■ ,x^} C D. The right hand sides in (Sal and (8b I become a vector 
and a symmetric positive semidefinite K-hy-K matrix, respectively. Hence a(x) becomes 
a vector of tabulated values for x € Dk- It is natural to have one Wiener process per 
point x'' in the spatial discretisation, so that K = M. The corresponding K tabulated 
individual diffusion coefficient functions, bj, will be obtained by a square root factorisation 
of the computed matrix, by means of an eigenvector expansion; this choice of factorisation 
preserves the connection between the evaluation point Xk and the elements k in bj and 
produces spatially localised functions, consistent with the association of individual Wiener 
processes and points in D^- 

In the second step, the initial configuration, Xq, in ^ is chosen so that the microscopic 
domain D includes a solid-liquid interface in equilibrium. Since the interface is stationary 
no phase transformation occurs in the simulation, and consequently the part of the reaction 



term in the AUen-Cahn equation ( 2b ) relating the speed of the phase change to the deviation 
from the melting point, k^k-^g' {<t>){Tyi — T), can not be obtained; the simulation must 
be performed at the melting point, Tm, under the given conditions. The simulation of 
a travelling front, off the equilibrium temperature, requires more advanced micro model 
simulations than the ones considered here. 

The interface is assumed to be locally planar on the microscopic scale and the spatially 
averaged properties are expected to vary much more slowly in the directions parallel to the 
interface than in the direction normal to the interface. Label the direction normal to the 
interface as direction xi and let X2, be orthogonal directions in the plane of the interface. 
Then the moUifier, r/, in ([5| can be chosen to make the averages much more localised in the 
Xi direction than in the X2 and a;3 directions. In the microscopic domain, D, the averages 
in the X2 and x^ directions are chosen to be uniform averages over the entire domain, so 
that the phase-fields, m and meg, and the drift and diffusion functions, a, /3j,fc, a, and bj, 
become functions of one space variable, xi. Hence the evaluation points in Dk are only 
distinguished by their xi coordinates. As mentioned above, the drift coefficient, a, depends 
on the derivatives up to second order of, 77, and the potential energy contributions mi. After 
averaging out the X2 and X3 dependence, it can be written as 

a(xi; X*) = k^T^m{xi-X') + —A^{xi-X') + ^0(2^1; ^*), 
for some functions Ai and Aq. Keeping this form in the averaging, the coarse-grained drift 



coefficient in (Sal can be written 

axi 0x1 



a{xi) = kBT^-^m^^{xi) + -^ai{xi) + ao{xi), 



where the second order derivative of the averaged phase- field. 



ma-v{xi) = 



T 

m(xi; X*) dt 







(9) 



corresponds to the diffusion term in ( 2b ) . Assuming that the averaged phase-field rriav is 
a monotone function of xi in the interface, the explicit dependence on the spatial variable 
can be eliminated by inverting mav and defining 

a(mcg) = a(TOa:v (meg)), 6j(mcg) = 6j(mav (meg)), (10) 
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which give drift and diffusion coefficients on the form ([7]). 

The present study is a practical test of the method described above, fn particular the aims 
are to verify that Smoluchowski dynamics can be used in practise, in the sense that the 
coarse grained drift and diffusion coefficient functions can be determined together with the 
phase-field model potential, /, and that they seem reasonable. For this purpose simulations 
are performed at just one temperature and density (at the melting point) and with just 
two values of the angle of the stationary interface with respect to the crystal structure in 
the solid. An actual determination of the model functions in the phase field model would 
require many more simulations with varying parameters. 



2 Computational Methods 

The numerical computations consist of molecular dynamics computations, giving the mi- 
croscopic description of the two-phase system, and the extraction of model functions for a 
coarse grained stochastic differential equation model. 



2.1 Molecular Dynamics Models and Simulation 

Two mathematical models of the material are used; both are one component molecular 
dynamics models where the interaction between particles is determined by a pair poten- 
tial of the exponential-6 (Exp-6) type. The coarse graining is based on a stochastic model 
where the particle trajectories on the diffusion time scale are given by the Smoluchowski 
dynamics The computations with this model are performed under constant volume at 
the melting point where a liquid and a solid phase coexist in the computational domain. 
The melting point is determined using constant pressure simulations of the deterministic 
molecular dynamics model where the particle trajectories are determined by Newton's sec- 
ond law with forces given the by gradients of the model potential. Both models and the 
corresponding simulations are described below, after a description of the potential common 
to the models. 



2.1.1 Pair Potential Defining the Total Potential Energy 

The microscopic system consists of TV identical particles at positions X — {Xi, . . . , Xj^) in 
three dimensions. The total potential energy, U, of the system is determined by the particle 
positions through 



using pairwise interactions only. The pair potential is the spherically symmetric Exp-6 
potential 




N 



N 



(11) 



<i)(r) = Aexp{-Br) 




(12) 
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with r denoting the distance between two particles, and A, B, and C being positive 
model parameters. The Exp-6 potential, like the similar Lennard- Jones pair potential, 
^Lj('') = 4eLj ((o'Lj/?')^^ — (ctlj/?")^), is a short range interaction that can be used to 
model condensed noble gases. With the parameters used here, obtained from [IT], the 
Exp-6 potential models Argon at high pressures. At pressures around 2 GPa, where the 
solid-liquid phase transition will be simulated, the Exp-6 potential with its slightly softer 
repulsive part describes the equation of state of Argon better than the Lennard-Jones po- 
tential does; see [TTJ [T5] . The shapes of the two pair potentials around the global minimum 
of the Lennard-Jones potential can be compared in Figure 2(a) the typical inter atomic 



distances between nearest neighbours in both the simulated solid and liquid will be close to 
1. Note that, while the Lennard-Jones pair potential tends to infinity as the interatomic 



Pair potentials 0(r) in reduced Lennard-Jones units 



Tile Lennard-Jones Pair potentiai, O (r) 



energy g 




(a) Pair Potentials for Argon 



(b) General form of the Lennard-Jones poten- 
tial 



Figure 2: (a): The Exp-6 pair potential is similar to the Lennard-Jones pair potential near 
the minimum, but the repulsion is slightly weaker in the Exp-6. The radius and the energy 
are measured in reduced Lennard-Jones units, where the Lennard-Jones parameters are 
CLj = fee 120 K and ctlj = 3.405 A. 

(b): The parameter ctlj is the radius where the Lennard-Jones potential is 0, which is 
equal to the potential at infinite separation, and the parameter elj is the depth of potential 
minimum. 



distance tends to zero, the Exp-6 pair potential, as stated in ( 12 1, reaches a global maximum 
before turning down and approaching minus infinity in the limit. This clearly illustrates 
that the model based on the Exp-6 potential breaks down if two atoms come too close, 
but neither one of the pair potentials is designed to describe interactions of particles much 
closer than the typical nearest neighbour separation. 



For short range potentials, like the Exp-6 and the Lennard-Jones potentials, the potential 
(and its derivative) decay sufficiently fast for the combined effect on the total potential 
energy (and the interatomic forces) of all atom pairs separated more than a certain distance 
to be negligible compared to the effect of the pairs separated less than the same distance. 
To take advantage of this in computations a cut-off radius is introduced and all interactions 
between particles separated by a distance larger than the cut-off are neglected; instead of 
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summing over all k^i in the inner sum in (11) the sum is only taken over particles in a 
spherical neighbourhood of particle i. 



All the physical quantities in this report are given in the reduced Lennard- Jones units. Thus 
length is measured in units of ctlj, energy in units of elJj and time in units of yj maf^j/ei^j, 
where m is the mass of one atom. (The time unit is the inverse of the characteristic 
frequency.) A list of the dimensionless units in the Argon model as well as the parameters 
in the Exp-6 potential can be found in Table 2.1.1 At the temperatures and pressures 



considered here, the stable phase of the Exp-6 potential is either the Face Centered Cubic 
(FCC) lattice or a liquid phase. 



Quantity 


Unit 




Constant 


Value 


Energy 


1.6568 • 10-2^ J 




/cb 


1 


Time 


2.1557-10-12 s 






Mass 


6.6412 • 10-26 kg 




Parameter 


Value 


Length 


3.405 • 10-1° m 




A 


3.84661 • 10^ 


Temperature 


120 K 




B 


11.4974 


Pressure 


4.1968 • 10^ Pa 




C 


3.9445 



Table 1: Atomic units and corresponding values of physical constants and parameters in 
the Exp-6 model (12). Non dimensional molecular dynamics equations are obtained after 
normalising with the atom mass, m, and the Lennard- Jones parameters, ctlj and eLj; in 
this Argon model m = 6.6412 • 10-26 kg (or 39.948 atomic mass units), ctlj = 3.405 A, and 
elj/^b — 120 K, where k-Q is the Boltzmann constant. 



2.1.2 Newtonian System Simulated at Constant Pressure 



The purpose here is to approximately determine the melting point at a high fixed pressure, 
to be able to set up and simulate stationary (FCC-liquid) two-phase systems later. Deter- 
mination of the melting point follows the two-phase method described by Belonoshko and 
co-authors in pj. 



The mathematical model is a classical system of TV identical particles where the po- 



sitions, X* — {X{, . . . , X^f), and the velocities, z)* 
to Newton's equations 



{v\, . . . , w^), evolve in time according 



~dt~ 
~dt 



(13a) 
(13b) 



where the total potential energy of the system is given by (11)-(12| using the parameter 
values in Table 2.1.1 Here V x denotes the gradient with respect to the particle positions. 
The force acting on particle i is —Vx U{X''^) and, since all particles have unit mass in 
the non-dimensional units, the acceleration is equal to the force. Particle positions are 
restricted to a finite computational box with periodic boundary conditions, corresponding 
to an infinite system where the same configuration of particles is repeated periodically in all 
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three directions; a particle leaving the computational cell on one side enters the cell again 
from the opposite side and particles interact with periodic images of particles in the cell. 



For a fixed volume of the computational cell the equations ( 13 1 will preserve the total 
energy, E, (the sum of potential and kinetic energy) of the system as well as the number 
of particles. It will approximately sample the {N, V, E) ensemble. In the determination 
of the melting point the simulations are instead performed in an approximation of the 
(A'', T, P) ensemble, using a constant number of particles, N, a constant temperature, T, 
and a constant pressure, P. This must allow for the volume of the computational cell to 
change during the simulation. There must also be mechanisms for keeping the temperature 
and the pressure constant, thus modifying (13 1 so that the total energy varies. 



Numerical computations of the (N,T,P) molecular dynamic simulations were per- 
formed using Keith Refson's publicly available software package Moldy, [9\. Constant 
temperature was enforced using the Nose-Hoover thermostat, where the equations of mo- 



tions ( 13 1 are modified, and extended, to include an additional degree of freedom mod- 
elling a thermal reservoir. The fictitious inertia associated with the thermal reservoir was 
100 kJmoP^ps^, corresponding to 21.57 in the dimensionless equation. The pressure was 
kept constant using the Parinello-Rahman equation, controlling the dynamics of the vec- 
tors (three edges) that define the computational cell. The fictitious mass parameter in the 
Parinello-Rahman equation was 300 amu corresponding to 1.20- 10"* in the reduced Lennard- 
Jones units. A short description of the Nosc-Hoovcr thermostat and the Parinello-Rahman 
equation, with references to papers with theoretical foundations of the methods, can be 
found in the manual [10 . 

The time stepping method in Moldy is a modification of Beeman's algorithm using predictor- 
corrector iterations in the computation of the velocities; see (TUj for details. The simulations 
described here used the constant time step 4.639 • 10~^ and the potential cut-off 2.937. 



In the two-phase method for determination of the melting point the molecular dynam- 
ics simulation starts from an initial configuration that is part solid and part liquid. As the 
{N, T, P) simulation proceeds the whole liquid part will solidify, if T < Tm for the given 
pressure, or the solid will melt, if T > Tm, resulting in a single phase system. Starting 
from a coarse estimate of the temperature interval containing the melting temperature, 
that interval can be narrowed down by running simulations at temperatures in the interval 
and noting whether they equilibrate to an all solid or an all liquid system. The validity of 
this two-phase approach has been verified in [1 for determining, among other things, the 
melting point of a molecular dynamics model of Xenon, similar to the Argon model used 
here. 

The initial configuration in a two-phase simulation was composed of pre-simulated solid and 
liquid configurations. The solid part was prepared by taking a perfect FCC configuration 
and performing a short molecular dynamics run at the temperature and pressure of the in- 
tended two-phase simulation to adapt the size of the computational cell. Initially the sides 
of the computational cell were aligned with the sides of the unit cube in the perfect FCC lat- 
tice; see Figure [3] While in general the dynamics of the cell edges in the Parinello-Rahman 
equations allow the cell to take the shape of any parallelepiped, here the dynamics were 



10 



(a) The FCC unit cube 



(b) Eight FCC unit cubes 



Figure 3: A perfect FCC lattice consists of FCC unit cubes, (a), stacked next to each other 
in three dimensions, (b). With one atom in the (0,0,0) corner of the unit cube the three 
other atoms are placed at the centres of the cubic faces intersecting in (0,0,0). 



restricted to only allow rescaling, without rotation, of the three edges and thus keeping the 
rectangular box shape of the cell. The preparation of the liquid part started from the con- 
figuration of the already prepared FCC-solid and a run was performed at a temperature well 
over the estimated melting point, where the sample would melt quickly; after equilibrating 
at the higher temperature the sample was quenched to the temperature of the two-phase 
simulation. Only one side of the computational cell was allowed to change while preparing 
the liquid part and thus the orthogonal cross section of the simulation cell was preserved 
from the FCC simulation. The solid and liquid parts were joined in the two-phase initial 
configuration by placing them next to each other, letting the cell faces of identical shape 
face each other. The general appearance is similar to the configurations shown in Figure [5] 



on page 16 even though those configurations belong to the constant volume Smoluchowski 



simulations where the set up procedure is slightly modified. Periodic boundary conditions 
were still applied in all directions, so that each part (solid or liquid) corresponded to a 
semi-infinite slab surrounded on two sides by the other phase with the effect of simulating 
a periodic, sandwiched, material. Voids of thickness of approximately one nearest neigh- 
bour separation were introduced in both solid-liquid interfaces to make sure that no pair 
of particles ended up to close in the initial configuration. Since the two-phase simulations 
were performed at constant pressure, the voids would fill in the beginning of the run as the 
length of the computational cell decreased. 



In the two-phase simulations the lengths of all three vectors defining the cell edges were 
allowed to change. Starting from an initial two-phase configuration the molecular dynam- 
ics simulation was run until the system was considered equilibrated. After equilibration 
the computational cell was filled with either the solid or the liquid phase. The density of 
the FCC solid is higher than that of the liquid phase. If the phase change was solidifi- 
cation of the liquid, then the volume of the computational cell would decrease during the 
equilibration stage before assuming an approximately constant value; if the solid was melt- 
ing, the total volume would grow during equilibration. The density of the stable phase at 
the given pressure and temperature was obtained by time averages of the simulation after 
equilibration. 

When the volume per particle is shown as a function of the temperature, at constant 



11 



pressure, it will display a sharp change at the melting point; see Figure 4(a) on page 12 The 
procedure will obtain an interval around the melting point and the accuracy can be improved 
by performing simulations at more temperatures to shorten the interval of uncertainty. 
However, the equilibration requires longer time when close to the melting point and the 
cost for refining the approximation grows, not only because the number of simulations 
grows, but more importantly because every single simulation takes longer to perform. 



Simulation at pressure 47.6554 



E 
o 

ra 0.8 



o 

> 




Two pfiase system at 
T = 2.9 and V = 0.7883 



0.76 
2.5 



2.8 2.9 3 3.1 
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(b) Extrapolated data at T = 2.9 



Figure 4: (a) Volume per atom as a function of temperature at a the pressure 47.6554 
(or 2.0 GPa) in (iV, P, T) simulations. Data points from simulations that are considered 
equilibrated are marked with o and those from simulations that are not equilibrated are 
marked with x. The two regions of equilibrated values where the volume per atom varies 
approximately linearly correspond to solid (FCC), at lower temperatures, and liquid, at 
higher temperature, respectively. The melting point at the given pressure is somewhere 
in between; the approximate value T — 2.9 is used below and in the constant volume 
simulations. 

(b) The volume per atom of solid and liquid have been extrapolated to T — 2.9 by least 
square fits of straight hues to the simulation data and the corresponding number densities, 
p, have been computed. If T = 2.9 is sufficiently close to the melting point at this pressure, 
then the two phases will coexist in constant volume, {N,V,T), simulations provided that 
the total density is between the estimated densities of pure solid and pure liquid. The 
ratio of the volumes of the solid and the liquid part is determined by the total density 
of the combined system. The tabulated value of the density for a combined system gives 
approximately equal volumes of both parts at a pressure close to the one in the constant 
pressure simulations. 



The main purpose here is to investigate the possibility of obtaining the model functions 
in a coarse grained phase-field model from {N,V,T) Smoluchowski dynamics simulations. 
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as described next. Therefor the accuracy in the determination of the melting point at the 
given pressure is critical only to the extent that it must be possible to perform the constant 
volume simulations at this temperature; that is, it must be possible to perform simulations 
on a two-phase system with stable interfaces between the solid and liquid parts. If the 
purpose were to perform computations at the melting point at this very pressure, then 
more computational effort would have to be spent on the accuracy of the melting point and 
the corresponding densities. 



The numerical simulations were performed with N — 8000 particles; the initial solid con- 
figuration consisted of 4000 particles, corresponding to 10 x lOx 10 FCC unit cells with four 
atoms each, and the liquid had the same number of particles. From simulations at the 
pressure 47.7 in the reduced Lennard- Jones units (corresponding to 2.0 GPa) an approx- 
imate value of 2.9 for the melting point was obtained together with number densities for 
the liquid and solid extrapolated to this temperature; see Figure |4] on page 12 Fixing the 
temperature and the number density N/V, only one degree of freedom remains in the triple 
{N, V, T) , allowing the system size to vary. 



2.1.3 Smoluchowski System Simulated at Constant Volume 



The constant volume and temperature Smoluchowski dynamics two-phase simulations de- 
scribed here were used to compute the functions (10) defining the coarse-grained phase- 
field dynamics ([t]), as described in the introduction. This meant computing time averaged 
quantities like the time averaged potential energy phase-field ^ and the corresponding 
coarse-grained drift and diffusion coefficient functions ([s]). 



The mathematical model is that of N particles whose positions X* follow the Smolu- 
chowski dynamics 

dX* = ~\7xU{X^) dt+y^2kBTdW\ (14) 

introduced on page [3] There are no velocities in the Smoluchowski dynamics. Instead 
the positions of all particles in the system give a complete description of the system at a 
particular time. Such a description, X*, will be refered to as a configuration of the system. 
The particles are contained in a computational cell, shaped like a rectangular box, of fixed 
dimensions and the boundary conditions are periodic in all directions. Hence the volume, 
V, and the number of particles, N, are fixed. Without velocities there is no kinetic energy, 
but the temperature, T, enters directly in the dynamics. The temperature parameter is 
held fixed, which can be viewed as a kind of thermostat built into the dynamics. 

Since the volume of the computational cell is constant, unlike in the {N,T,P) simulations 
above, the overall density of the system remains constant over time, which allows for sta- 
tionary two-phase configurations where part of the domain is solid and part is liquid. 



The numerical simulations The discrete time approximations X" of X*", were com- 
puted using the explicit Euler-Maruyama scheme 

X" = - Vx C/(X"-i) At" + ^/2kBT AW"-, (15) 
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where At" — i" — is a time increment and AVF" — — is an increment in 

the 3iV-dimensional Wiener process. Each run was performed using constant time step size, 
Ai" = At, but the time step could change between different runs depending on the purpose; 
in the equihbration phase the typical step size was At = 10^*, but in the production phase 
the step size had to be taken smaller, as discussed later. 

The computation of VxC^(-^"~^) in every time step is potentially an 0{N'^) operation 
since the potential is defined by pairwise interactions. The computations described here 
used the potential cut-off radius 3.0, which meant that each particle only interacted directly 
with a relatively small number of neighbours (independent of N since the density was 
approximately constant). To avoid the 0{N^) task of computing all pairwise distances 
in each time step, the computational cell is divided into smaller sub cells, where the size 
is defined in terms of the cut-off radius so that two particles only can interact if they 
are in the same sub cell or in two neighbouring sub cells; information about particles 
migrating between sub cells is exchanged in each time step. The computations use a two 
dimensional grid of sub cells, where the particle positions within each sub cell are sorted 
with respect to the third coordinate dimension in every time step. When the particles are 
sorted the sweep over all particles in a sub cell can be efficiently implemented and the sorting 
procedure is not too expensive since the particles do not move far in one time step. A more 
thorough description of this algorithm can be found in [13] . The actual code used here is a 
modification of a parallelised code for Newtonian molecular dynamics obtained from Mans 
Elenius in Dzugutov's group [1]; the main modifications when adapting to Smoluchowski 
dynamics is the removal of velocities from the system and the introduction of a pseudo 
random number generator for the Brownian increments, AM^". 

With the cut-off radius 3.0 used in the computation and the model parameters in Table [2TO] 
on page [9] the Exp-6 pair potential and its derivatives are small at the cut-off radius. Still 
the potential will be discontinuous at the cut-off, unless it is slightly modified. A small 
linear term is added to make the potential continuously differentiable at the cut-off radius. 
In the practical computations, both the pair potential and the derivatives were obtained by 
linear interpolation from tabulated values. 

The random number generator for normally distributed random variables was the Ziggurat 
method, described in [5], in a Fortran 90 implementation by Alan Miller, accessible from 
Netlib jj. The underlying 32-bit integer pseudo random number generator is the 3-shift 
register SHR3. Since the purpose of the simulations only is to investigate if the coarse- 
graining procedure gives reasonable results just one pseudo random number generator was 
used, while several different random number generators ought to be used in a practical 
application. The generator was initialised with different seeds on different processors in the 
parallel computations, but it does not have distinct cycles simulating independent random 
variables. The hope is that the nature of the molecular dynamics simulations is enough 
to avoid the danger of correlated random numbers on the different processors, but this 
could be tested by comparing with other pseudo random generators that actually simulate 
independent random variables on different processors. 



The two-phase systems for the Smoluchowski dynamics simulations were set up to 
obtain a two-phase system at temperature T — 2.90 with approximately equal volumes of 
solid and liquid and with stationary interfaces. To achieve this two equal volumes of FCC- 
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solid and liquid were pre-simulated with the densities tabulated in Figure |4] on page [T2j 
The preparation of the initial configurations for the Smoluchowski dynamics two-phase 
simulations was similar to the procedure described above, but some adjustments must be 
made because of the constant volume restriction. The shape of the computational cell used 
when generating the solid part was chosen to match the periodic structure of the FCC 
lattice at the tabulated density for the FCC part. A short equilibration run, at T = 2.90, 
starting from a perfect FCC lattice at this density gave the initial solid configuration. The 
computational cell for the initial liquid part was chosen to be the same as the one in FCC 
simulation and the initial configuration when pre-simulating the liquid part was obtained 
from the FCC configuration by distributing vacancies to get the correct density in the liquid. 



In a simulation of (151 using a temperature, T, above the melting point, Tm, the sample 
was melted and equilibrated. Afterwards the liquid was cooled to desired temperature using 
a subsequent simulation with T = Tm ■ 

Since no pair of atoms can be too close in the initial configuration, gaps had to be intro- 
duced between the solid and liquid parts, but the voids could not be introduced as additional 
volumes in the computational cell; the individual parts were equilibrated at (iV, V, T) corre- 
sponding to the expected densities for solid and liquid in the combined system, so increasing 
the total volume would reduce the overall density, resulting in partial or total melting of the 
solid part. To make room for the voids both the solid and the liquid parts were compressed 
slightly in the direction normal to the solid-liquid interfaces, before inserting them in their 
respective volumes in the computational cell for the two-phase simulation. Initial configu- 
rations obtained by this procedure are shown as configurations (a) and (c) in Figure [s] on 



page 16 The orientation of the solid-liquid interfaces with respect to the FCC lattice differ 
between the two initial configurations shown, and these orientations with the corresponding 
numerical simulations will be labelled Orientation 1 (01) and Orientation 2 (02) in the 
following. The shaded plane in Figure [6|^b) shows the orientation of the interface in 01 and 
the shaded plane in Figure ^c) shows the orientation in 02. 

Even though the compression in one direction was small, it introduced an artificial internal 
stress in the system. The higher value of the phase-field in the subfigures (a) and (c) in 
Figure [s] compared to the corresponding regions in the subfigures (b) and (d) is an effect 
of the compression. In the initial phase of the equilibration of the two-phase system, the 
compressed parts expand to fill the voids. The phase-fields in the interiors of the solid and 
liquid parts in subfigures (b) and (d) have reached the levels seen in the corresponding 
single phase systems, which shows at least that the local potential energy contributions had 
returned to normal before the production runs started. 

As a test of the two-phase configuration serving as initial data in the production run, the 
radial distribution functions in the interior of the two phases were computed. The radial 
distribution function, g{r), is useful for identifying the phase of a single-phase system. For 
a single component system g{r), where r £ R"*", is implicitly defined by the condition that 
the average number of atoms in a spherical shell between the radii ri and r2 from the centre 
of any atom is 

/•r2 

P / .g(r)47rr^ dr, 
J n 

where p is the global particle density. In other words, the radial distribution function is the 
average particle density, as a function of the separation r, normalised by overall density. 
Figure [t] on page 18 shows good agreement for simulation 02 between g{r) corresponding 
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(a) Initial configuration, orientation 1 



1 ,0 '- 




(d) Configuration at a later time, orientation 2 



Figure 5: Snapshots of the process of setting up initial configurations for the two- phase 
simulations 01 and 02. The left part is solid (FCC) and the right part liquid. In the 
initial configurations, (a) and (c), the individual parts have been equilibrated at T^^it (for 
the combined system), and slightly compressed in one direction (to allow for two gaps). 
Subfigures (b) and (d) show configurations at later times when the parts have expanded 
to fill the voids and form two interfaces. The atoms are coloured according to a computed 
phase variable; in (a) and (b) the phase variable is just the instantaneous field m{xi; X'^), 
whereas (b) and (d) use discrete time averages approximating J^^ m{xi; X*) dt. 

Simulation 01 used 64131 particles in a computational cell of dimensions 
93.17 X 23.29 x 23.29, while simulation 02 used 78911 particles in a cell of dimen- 
sions 100.86 X 24.71 X 24.96. 
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(a) Unit cell 



(b) Orientation 1 



(c) Orientation 2 



Figure 6: The shaded planes in (b) and (c) show the two orientations of the sohd-hquid 
interface with respect to the FCC lattice treated in the numerical simulations. 



to single phase solid and liquid configurations and g{r) computed in the interior of the two 
phases, excluding two intervals of length 10.0 in the interface regions. 

An effect of the finite size of the computational cell is that periodic boundary conditions 
may interact with the solid and affect the results; here the computational cell was chosen to 
match the FCC structure in a specific orientation with respect to the box and thus stabilises 
the structure and orientation. It is important to know that the density in the FCC part (and 
hence the box cross section) is consistent with constant pressure simulations close to the 
melting point. A related question is whether the length of the computational box is large 
enough for properties around the interfaces in the infinitely layered structure to be good 
approximations of those near an interface between a solid and liquid on the macroscopic 
scale. 



2.2 Computation Of the Coarse-Grained Model Functions 



The coefficient functions ( 10 1 in the stochastic differential equation ([7| for the coarse-grained 
phase-field are defined in terms of the time averaged expected values ([s]) and ^ on the 
form 











Jq 





where is a configuration of a stationary two-phase system. By setting up an initial 
configuration, Xq, as described in the previous section, and simulating discrete sample 
trajectories using the Euler-Maruyama method (15 1, a sequence of configurations {X''}^ 



approximating the sequence {AT**}^]^ for some times < ti < • • • < = T, is obtained. 
In a post processing step a set of configurations S C {A'''}^^ is selected and averages 

xes 

consistently weighted with weights wx , are computed as approximations of the correspond- 
ing expected values in the continuous time model. It is usually more efficient not to include 
every configuration in the averages. This will be discussed in Section [3j 

As described in the introduction, the averages are functions of the coordinate direction xi, 
normal to the planar interface, since the moUifier in the definition ^ of the microscale 
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Radial Distribution Function, g(r) 



Pure solid phase 

Solid in two phase system 




Radiai Distribution Function, g(r) 

Pure liquid phase 

— Liquid in two phase system 
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(a) FCC 



(b) Liquid 



Figure 7: The radial distribution function, g{r), computed from several configurations, 
separated in time, in the process of setting up the two-phase system in simulation 02. The 
solid curve shows g{r) computed as an average over all particles in the computational cell 
used while pre-simulating the solid and the liquid part, in subfigure (a) and (b) respectively. 
The dashed curves show g{r) computed as an average over particles in two slices of the 
computational cell of the two-phase system; subfigure (a) shows g{r) obtained from the 
slice 5.0 < xi < 45.43, inside the solid phase, and subfigure (b) shows g{r) from the 
slice 55.43 < xi < 95.86, inside the liquid phase. The configurations are taken from an 
equilibration run, after the closing of the initial gaps between the pre-simulated phases, 
but before the "production" run. The radial distribution functions show good agreement 
between the single phase systems and the corresponding solid and liquid subdomains away 
from the interface. 
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phase-field, m, is chosen to take uniform averages in the planes parallel to the interface. 
The moUifier used in the computations is 



r]{x) = f]{xi) = cexp ^-^ (^^^ 



<-Rci 



(16) 



where c is a normalising constant, e is a smoothing parameter, and Rc is a cut-off. The 
smoothing parameter is on the order of typical nearest neighbour distances, e « 1, and 
Rc = 6e, for all choices of e, which gives r]{Rc) « 1.5 • 10~^77(0); the shape of -q can be seen 



in Figure 25(a) on page 41 



An explicit derivation of expressions for the drift and the diffusion is given in Appendix |Aj 
Separating the drift in terms containing two, one, and zero, derivatives of the moUifier, the 
right hand side of (Sal is approximated by 



d_ 
dxi' 



knTg^As (m) + ^^^^^5 (oi) + ^5 (ao), 



where 



N 



and 



a,ix;X)^J2(kBT^m,{X))[F,iX)]Mx~X^) 



^ / 1 \ 

ao{x;X) = ^bTVx, ■ F,{X) + -\\F,{X)\\^ ~ X^ 

j=i V / 

N N 

E h{X)-F,{XHx^X^). 



(17) 



(18) 



Here Fj is the total force acting on particle j, [Fj{X)]i is the xi-component of the force, 
and fij are the contributions from individual pairs. 



N 



F,ix) = -vx,uix)= ^ ci>'(||x,-x;i^ 



N 



II 



The right hand side in equation (8b I, for the coarse grained diffusion, is approximated by 

N 



where 

P]{x,y;X 



B{;-) = As I 2k^TY,{Pj{-r;X) + q,{;--X)) |, 



(19) 



e 

m,{X) 
m,{X) 



N 



[x~X^],ri{x~xM[F,{X)],T^{y-X^)+ ^ [U,{X)],'n{y - X, 

{v-x^,^^{y-x^)[{F,{x)\^^{x-x.;)+Y. \!^Ax)V^^^-x^ 
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and 

1 / ^ 
q,{x,y;X)^-iF,{XMx-X^)+ ^ f,,{X)7^ix ~ X,) 

N 

F,{X)rj{y^X^)+ ^ l,,{X)^{y - X.^) 



The functions As {4') £^re computed in a discrete set of points Dk — {x\}fLi along the xi 
axis of the molecular dynamics domain. This makes the computed components, ^5 (w), 
^5(01), and ^5(00), of the drift coefficient function vectors and the computed B a 
K-hy-K matrix. The individual diffusion coefficient functions bj are obtained by taking 

the square root of the computed diffusion matrix, B = B^^^ (B^^^)^ , and letting the j:th 
1/2 - . 

column of B define bj . While an exact computation would produce a symmetric positive 
semi definite matrix B, finite precision effects make some computed eigenvalues negative, 
but small in absolute value. In an eigenvector factorisation of B, let A denote a diagonal 
matrix with all eigenvalues of B and A_|_ a smaller diagonal matrix containing the dominant, 
possibly all, of the positive eigenvalues but no negative ones. Let V and V+ be the matrices 
of the corresponding eigenvectors. Then the square root of the matrix A_)_ is a real diagonal 
matrix which can be used in the approximation 

B = VAV'^ « y+A+y+^ = ft/+A+i/V+^) (v+A+^/^V+^Y =: BB^. (20) 



With one Wiener process Wj in the coarse-grained stochastic differential equation ([t]) per 
evaluation point, K = M, the component vectors, bj, of the diffusion in coarse-grained 
equation can be defined as the column vectors of the matrix B, to obtain 



M 

Yb.bJ « B. 



If two grid points, xi and yi, are further apart than twice the sum of the cut-off in the 
potential and the cut-off in the moUifier, then pj{x,y;-) and qj{x,y;-) is zero; hence a 
natural ordering x\ < x\ < ■ ■ ■ < xf of the grid points makes B a band matrix. The 



definition of _B in ( 20 1 preserves the connection between grid points and diffusion functions 



and the dominating terms in a tabulated vector bj are those of nearby grid points. 



3 Results 



This section describes results from numerical simulations performed to compute the coarse- 
grained model functions. The value of the smoothing parameter e in the moUifier is 1.0, 
unless another value is specified. 
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3.1 The averaged phase-field ~ As {m) 

The first observation is that during the time intervals of the molecular dynamics simula- 
tions, the interfaces between the solid and the liquid subdomains were sufficiently stable 
for the averaged potential energy phase-fields, .45 (m), to appear qualitatively right. The 
phase- field appears to have two distinct equilibrium values, corresponding to the solid and 
liquid subdomains, and the transitions between the two regions are smooth and occur over 
distances of a few nearest neighbour distances; see Figurejs] Figure [9(b)] shows that the com- 
putational cells in the molecular dynamics simulations are large enough for the phase-field 
in the interior of the two phases to attain values similar to the values in the corresponding 
single phase simulations. In simulations with a cubic, 23.29 x 23.29 x 23.29, computa- 
tional cell the gap between the phase-field levels in the solid and the liquid was significantly 
smaller, which indicates that the length of the computational cell can not be taken much 
smaller than in simulations 01 and 02. It is still possible that further increasing the size 
of the computational cell may affect the results. 



3.2 The averaged drift a ~ As 



a 



The average As {m) approximates the expected time average ^ . The next expected value 
to study is the one defining the coarse grained drift in (8a). In a stationary situation, where 



the interfaces do not move during the simulation and the averaged phase-field converges to a 
stationary profile, the average total drift in the stochastic differential equation describing the 
phase-field variable must converge to zero. Still the time averaged total drift corresponding 
to the simulation 02, whose averaged phase-field was discussed above, is far from zero; see 



Figure 10 The computed time averaged drift 



^5(a(x;X")) « -E 



a{x;X*)dt 



X'' = X. 



depends both on the length of the time interval where the average is computed, the number 
of configurations used in the average, and on the discrete approximation X" of X*"; these 
potential error sources must be analysed to explain the result. 



3.2.1 The effect of discrete time dynamics 



First consider the error associated with the discrete dynamics. The explicit form of the 
drift is derived for the continuous time mathematical model with the Smoluchowski dy- 
namics (14 1, and not the discrete time Euler-Maruyama dynamics (15 1 that is used in the 
numerical simulations. For a fixed size of the time step this means that, even if the state 
of the numerical simulation is stationary on the time scale of the simulation so that time 
averaged phase-field converges to an equilibrium profile, the time averaged total drift will 
not go zero because of the time discretisation error. Figure [TT] shows that the computed 
radial distribution functions, here from single phase solid configurations, are close when the 
time steps used vary from 10"^ to 10"''; still the larger time steps give average computed 
drifts As {a{x; X")) that are inconsistent with the observed time evolution of the average 
phase-field As (m(a;;X")). As shown in Figure 13 



the time step At = 1 • 10 gives an 
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(c) Orientation 2 (d) Orientation 2 




(e) Orientation 2 (f) Orientation 2 

Figure 8: Subfigures (a) and (c) show the potential energy phase-field, As (m), computed 
from simulations 01 and 02, respectively. Subfigures (b) and (d) show the correspond- 
ing spatially averaged paxticle densities. Subfigures (e) and (f) show the pointwise sample 
variance associated with the averages in (c) and (d). The thick parts of the curves show 
the computed functions in molecular dynamics cell. The thinner parts show the periodic 
continuations across the boundaries of the cell, marked by circles. The averages in sim- 
ulation 01, and 02, were formed over 1721, and 1775, configurations separated in time 
by 5 • lO"**; so that the total time from first to last configuration was 0.860, and 0.8875, 
respectively. The high frequency fluctuations are small after averaging on this time scale, 
but larger fluctuations remain in both phases. This suggests that the two phase system is 
not yet equilibrated. Still the computed phase-flelds appear qualitatively correct. 
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(a) Density Phase-Field, 02 



(b) Potential Energy Phase-Field, 02 



Figure 9: The computational cell in the molecular dynamics simulations must be sufficiently 
large for the infinitely layered structure to resemble a system with a single solid-liquid 
interface on the macroscopic scale. In simulation 02 the total length of the computational 
cell was 100.86; subfigure (b) shows that this was sufficient for the averaged phase-field, 
As (m), to obtain values in the interior of each phase that are similar to the functions, 
marked by thick curves, obtained in the single phase configurations simulated during the 
setup of simulation 02. 
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Figure 10: The average total drift. As {a), based on the same 1775 configurations from 



simulation 02 as As (rn) in Figure 8(c) is still dominated by large oscillations 
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average drift that oscillates between -100 and -250, even when the computed phase-field 
As (m(x;X")) is approximately constant over times of the order 10. For this reason, the 
time step used in simulations 01 and 02, generating configurations for the computation of 
As {m{x] X")) and As {a{x; ^")), was At = 5-10^^ , while the time step used in the setup 
of the initial configurations often was a thousand times larger. With this small time step 
the fiuctuations in the computed average drift outweighs the deviation from the expected 
zero mean; see Figure |10| 

The choice of the time step size At — 5 • 10^^ was guided by a rough error estimate, taking 
into account the maximal absolute value of second order derivatives of the Smoluchowski 
drift — U (X*^) when the nearest neighbours don't come closer than approximately 0.8, as 
indicated by Figurc[TT] Then the time step was adjusted so that the slow convergence of the 
time averaged drift in terms of T and the number of configurations, X", was the dominating 
error source in the results. This over-killing of the time discretisation error in the molecular 
dynamics wastes computer power and could possibly be avoided by more accurate error 
estimates, allowing a matching of the different error contributions. Using a reasonable 
number of grid points, K, in the computation of the drift coefficient _fC-vectors and the 



diffusion K-hy-K matrix B, in (19 1, the computational cost for obtaining B in particular, 
far exceeds the cost of actually making a time step in the molecular dynamics simulation. 
Hence the additional cost of over-killing the time step error is not very significant, provided 
that not every configuration in the time stepping is included in the averages As (w). As (a). 



and B. In the averages shown in Figure pi and Figure 10 for example, the configurations 
were sampled at time intervals 5 • 10^**, corresponding to 1000 time steps in the molecular 
dynamics simulation. 

A further improvement may be to incorporate finite step-size effects in the expressions for 
the components of the drift. The higher order derivatives of the pair potential attain large 
values when two particles come closer than 1; see Figure [12] Hence the time step must 
be taken very small for Ito's formula to be a good approximation of the dynamics of the 
discrete system. Instead of a direct application of Ito's formula in the derivation of the drift 



and diffusion terms in ( 26 1 and ( 27 1 on page 44 one could include higher order terms in the 



expansion to improve the accuracy of the computed drift. 



3.2.2 Dependence on the length of the time averaging interval 



Next consider the dependence of the computed coarse-grained drift coefficient function on 
the length of the time interval T. Introducing the time averaged drift over a sample path 
as 



At=^ I a{--X')dt, 



the coarse-grained drift (Sal is a = E[v47-]. The rate of convergence of a, as T ^ cxd, in 



the continuous time mathematical model can be estimated by integration of the stochastic 
differential equation ([6| for the phase-field m. Integrating from to T gives 



r N 



m{-;X'^) -m(' 



/ a{-,X')dt+ j Y.ilM\X 

Jo Jo .1 



(21) 



j=i fe=i 
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il Distribution Functio 



Radial Distribution Function, g(r] 





(a) g{r), FCC, different At 



(b) Detail: first slope 



Figure 11: The radial distribution function, g{r), computed using four different step sizes 
in a single phase FCC simulation. The difference between the curves is small (a), even if 
the one obtained for At — 10^** differs visibly from the others in the first peak (b). In spite 
of the good approximation in the radial distribution function, the larger step sizes give very 
poor results in the computed dynamics of to. 



*(r) 



d<I>(r)/dr 



d'*(r)/clr' 



Figure 12: The absolute value of the Exp-6 potential and its derivatives grow very quickly 



with decreasing r, in the range with positive g{r) in Figure 11(b) The potential and its 
two first derivatives using the model parameters in Table |2.1.1[ on page[9l are shown here. 
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Drift: T cPm/dxJ +da^/dx^ +3^ 
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Figure 13: Using the step size At = 1 • 10""'' in the Euler-Maniyania scheme, the computed 
average phase-field ^5 (m) is approximately stationary during the time interval of the 
averaging. In subfigure (a) the average is based on 123 configurations, sampled at every ten 
thousandth time step, corresponding to a total time interval of 12.3. Still, the computed 
average drift As {ct) is far from zero during this time interval. The large deviation from 
zero is entirely due to the term As (ao)- 
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(a) Mean based on 111 configurations, T = 0.0555 (b) Mean based on 444 configurations, T = 0.2220 
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(c) Mean based on 1775 configurations, 7~ = 0.8875 
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1.9 (3.8 • 10^) 



(d) The spatial (xi) mean and, within parentheses, variance of ^5 (a) 



Figure 14: The total drift As {a), decays slightly faster with T than the predicted XjyT 
in the examples (a), (b), and (c) above. Here the number of configurations in the averages 
grows with T and the means and variances of As (a) tabulated in (d) suggest that the 
number of configurations still restricts the rate of convergence. The average in subfigure (c) 



is based on the same 1775 configurations from simulation 02 as As (jn) in Figure 8(c) 
The averages in subfigures (a) and (b) are based on the first 111 and 444 configurations 
respectively. 
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so that, by taking the expectation and using that, since X* is lF*-adapted, the expectations 
of the Ito-integrals vanish 



E 



r 



dt 



= E[m(.;X^)-m(.;X")] 



(22) 



Hence, if the phase-field is stationary, then the expected mean drift over time is zero. 



Normahsing (21 1 and (22 1 by T 



- E \At\ = 



( N 3 T 

- m{-,X^) - m(.; X°) - E [m{-,X^) - m(.; X°)] - ^ ^ / 



r 

and the variance of Aj- is obtained as 

Varpr] = E \{At - E [At])^ 
1 



(3,M-,x')dwl, 



T2 



Var 



m{-;X') -m(-;X") 



TV 3 



r2 



T2 



E 



i(.; X^) - m(.; X")) £ /3„.(-; ^*) dW^l^ j 



where last expression was simplified using the independence of the different components of 
W^, and the zero expected value of Ito integrals. Assuming that both the phase-field and 
all the diffusion coefficients are bounded, the dominating term in the expression for the 
variance is 



-I N 3 

^EEe 



T2 



i=l fc=l 



/3,-fe(-;X*) dWl,, 



In the two phase simulations considered here, the values of the computed phase-field varies 
between a lower level in the solid a higher in the liquid. Because of the small positive 
probability for two particles, with trajectories computed using the Euler-Maruyama dy- 



namics ( 15 ), to get within an arbitrarily small distance of each other, there is no guarantee 



that computed phase-field always will stay in this range. However, if the minimum inter- 
atomic distance becomes to small, that is a breakdown of the whole microscopic model 
and not just a problem when computing the drift; this situation has not been observed 
to happen in the simulations here and the observed values of the phase-field are all in the 
range (—1.5, 1.0). Hence the assumption that m is bounded seems reasonable here; a bound 
on the absolute value of the diffusion coefficients Pj ^ is less certain, and it will have to be 
larger than the bound on m. 



For the average drift to be small compared to the stationary values of the phase-field itself, it 
must be at least a factor 100 smaller than the computed average shown in Figure [lO] Based 
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on the rough analysis above, the expected time average of the total drift can be expected to 
decay as I/a/T with a large constant factor. When the computed drift As (a) in Figure 10 



is compared to averages computed using two smaller subsequences of configurations, the 

Even when 



14 



convergence to zero appears to be slightly faster than l/y/T; see Figure 
extrapolating with the measured convergence rate, decreasing the average drift by a factor 
100 would require increasing the averaging time interval by more than a factor 1000, which 
is beyond reach within the present project. With increasing accuracy in the time average, 
eventually the time step in the molecular dynamics simulations must be decreased, further 
increasing the computational cost. 

Since the total drift coefficient function, a(xi) ~ As {a{xi; •)), where 



d 



As (a(xi; •)) = k^T^-^As {m{xi; •)) + ^zT^s {ai{xi; ■)) + As (ao(xi; •)), (23) 



in the coarse grained model is expected to be zero in a stationary situation, a more accurate 
computation would serve primarily as a consistency test. On the other hand, the individual 
terms in the right hand side are not all expected to vanish independently. Indeed, it is clear 



from the results on As {^{xi] ■)) in Section 3.1 that the term with two differentiations with 



respect to xi will not be identically zero. This also shows that while the total drift is far 
from As (^(a^i; ■)) converged, at least one term is reasonably accurate. 

A closer look on the terms of the drift, reveals that the different terms are of different 



orders of magnitude. The term As (ao(a^i; •)); with defined in (18), contains both second 
order differentials of the potential with respect to the particle positions and second powers 
of first order differentials. These terms, as illustrated in Figure [12] attain much larger 
values than the potential itself and cancellation is required to reduce As {clq{xi] •)) to a 
size comparable with the two other terms in the drift. Figure |15(e) shows an individual 



ao(a;i;-) computed from one configuration; in the length of the computational cell, the 
values range from approximately -500 to +500, whereas the phase-field, m(xi; •), is of the 
order 1, and ai{xi; ■) is of intermediate magnitude. A comparison between the computed 



averages As {oi{xi; •)) in Figure 14 and As (ao(a;i; •)) in Figure 15 shows that As (ao(a;i; •)) 



is the dominates the other two terms completely here. 

The average As {o,i{xi; •)), contains first order differentials of the potential, but only to 
the first power. The convergence of is faster than that of As {<io{xi \ •)), but the computed 



averages in Figure 16 still show significant fiuctuations. The final term in As {oi{xi] •)) is 



kBT-^As {fn{xi; •)), which only depends on the potential and not its derivatives. This 
average converges faster than the other two and, even after two differentiations with respect 
to Xi, the fluctuations are small compared to the distinct structures at the interfaces; see 
Figure [17] 



3.2.3 Obtaining the phase-field double-well potential from the drift 



When defining a phase-field variable in terms the potential energy in the microscale model 
in Section JTl the goal was to compute a reaction-diffusion equation, like the Allen-Cahn 
equation (2b), for the coarse-grained phase-field. In a one dimensional problem, with T = 
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(a) Mean based on 111 configurations, T = 0.8875 (b) Mean based on 444 configurations, T = 0.8875 
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(c) Mean based on 444 configurations, T = 0.2220 (d) Mean based on 1775 configurations, T = 0.8875 

„ .... ^ Variance of drift component a„ 

Drift component a - instantaneous .-^ ^ a 
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(e) One single configuration, X* 



(f) Variance based on the configurations in (d) 



Figure 15: The term As (ao) is the slowest converging average in the drift average; a 
comparison with Figure [14] shows that this term dominates the total drift average. This 



explicit form of the term, given in (18) is a sum over all particles of terms that are second 



order in the particle forces and a term containing the divergence of the particle force; in 
the molecular dynamics simulation, these terms are large and so is the function oq, when 
computed from a single configuration, as in (e). Eventually the average must decrease to 
order 1 through cancellation, but for the number of configurations available here fluctuations 
dominate the computed averages As (ao)- 
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Drift component 
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(a) Mean based on 111 configurations, T = 0.8875 (b) Mean based on 444 configurations, T = 0.8875 



Drift component a^ 
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(c) Mean based on 444 configurations, T = 0.2220 (d) Mean based on 1775 configurations, T = 0.8875 



Figure 16: The term As (ai) is supposed to approach zero as the number of configurations, 
and T, increases, provided that the interfaces are stationary. Though the fiuctuations are 
large here, they are much smaller than in Figure [15] When the fiuctuations decrease a 
pattern appears with peaks at the two interfaces. This supports the observation, from the 
computed As {^n) in Figure [s] that the two phase system is not in equilibrium yet and the 
interfaces are not really stationary on the time scale of the average. 



31 



Drift component T cPm/dx^ Drift component T cf'm/dx^ 




X, X, 

(a) Mean based on 111 configurations, T = 0.8875 (b) Mean based on 444 configurations, T = 0.8875 



kgT d^m/dx^ Drift component k^T d^m/dx^ 
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(c) Mean based on 444 configurations, T = 0.2220 (d) Mean based on 1775 configurations, T = 0.8875 

Figure 17: The average k-oT^^As (cti) converges faster than the two other terms in As {ex). 

The fluctuations arc larger in subflgurc (c) than in (b), which indicates that the error is 
dominated by the length of the averaging time interval rather than the number of configu- 
rations sampled within the time interval. 
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Tm and ki constant, the AUen-Cahn equation reduces to 

dt ^ dx\ 

where the derivative of the double-well potential / gives the reaction part in this reaction 
diffusion equation. Now, the coarse-grained equation 



~ ^2/'(0) + noise, (24) 



dmlg{xi) = ( kBT-^mlg{xi) + —ai{xi) + ao(xi) j dt + y^bj{x) dW. 



M _ 

't 
j ' 



where 

ai{xi) ^ As {ai){xi), ao(a;i) = ^5 (ao)(a;i) JorxieDx, 



and the diffusion coefficient vectors, bj, are obtained from the factorisation (20 1, is a stochas- 
tic convection-reaction-diffusion equation. As the described above the time averaged drift 
is zero in a stationary situation, but in the computations presented here the fluctuations are 
still too large. In the ideal situation for a stationary interface, when all three components 
in the drift average have converged, the convection should vanish, that is 

OX I 

and the reaction and diffusion parts should cancel each other, so that 

= ^Br^m*g(xi)-Hao(a;i). (25) 

The second best thing, when some of the computed averages contain too large errors, is 
to extract information from the most accurate part, that is fceT J|^rnav(a;i). Assuming 
that this computed average already is close to what it would be in the ideal situation, an 
approximation of the reaction term can be obtained from (|25|) . 



The expression of the drift in the coarse-grained equation (|7]) as a function of the coarse- 
grained phase-field meg in the interface regions, instead of the space variable xi, assumes 



monotonicity of the phase-field near the interfaces to allow the inversion in (10). Figure 18 

dx 



shows mav(a;i) and fcB2^^^-^av(a^i) in the interval of monotonicity for mav(a;i) in the 



simulation 02. Using the computed fceT J^TOav(a^i) in (25 1, gives 



ao(xi) = ~kBT-^ms,v{xi). 



dx\ 

Inverting the computed function TOav(2^i) in the interface intervals, the derivative of the 
double-well potential / can be identified as 

/'(meg) = ao(mav (meg))- 
Integration with respect to meg in the interval between meg, ,• . and meg,- • . gives the 



double-well potentials shown in Figure 19(a) As expected the potentials obtained from the 



two different simulations 01 and 02 are slightly different. However, the potentials obtained 
from the two different interfaces in one molecular dynamics simulation cell also differ slightly 
and it is not possible to say that difference between simulations 01 and 02 depend on the 
orientation of the interfaces with respect to the crystal lattice. The computed double wells 
seem to be qualitatively right. 
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0.6 r 




-10 10 20 30 40 50 60 

Figure 18: The computed rrig^v^xi) in its monotone intervals in the interfaces together with 
the corresponding diffusion part of the drift kQT-^As (o-i)- The curves shown are part of 
the those in Figure [8(d)] and Figure [l7(d)| 




Figure 19: 

(a) The computed double well potentials from both simulation 01 and 02 using TOav shown 
in Figure jsj and the corresponding A;BrJ^mav(a;i). 

(b) The computed double well potentials from one of the interfaces in 02, using three 
different values of the smoothing parameter e in the moUifier. Since the interface width 
varies with e the height of the potential barriers vary with e. Here double- wells have been 
rescaled with factors obtained in the analysis of the e-dependence in Figure 27 to compare 
the shape of the curves. 
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3.3 The averaged diffusion matrix B and the coarse-grained diffu- 
sion coefficients b^. 



The final component to extract in the coarse-grained model is the diffusion in the stochastic 
differential equation for m*^. Using e = 1.0 and the same 1775 configurations that were used 
in the computation of the averaged phase-field and drift for simulati on 02, the averaged 
diffusion matr ix B, has been computed, with the result shown in Figure |20(a)[ As described 



in Section 2.2 the square root of Bis computed by an eigenvector decomposition where all 
negative eigenvalues are set to zero; the result is shown in Figure |20(b)[ The negative 
eigenvalues are very small in absolute value, compared to the dominating positive ones, 
so the error made by neglecting them is insignificant when BB^ is compared to B. By 
choosing the diffusion coefficients bj in the coarse-grained stochastic differential equation 
as the columns of B, they become localised in space; see Figure 20(c)| With e = 1.0 the 
observed difference between the diffusion in the solid part and the liquid part is small, as 
shown in Figure |21[ 



3.4 Dependence on the smoothing parameter 



The moUifier r/ includes a parameter, e, determining the scale on which the local average is 
taken. This is in itself an ad hoc variable in the micro model and it is important to analyse 
its effects on the computed quantities. 



A lower limit on e is set by the demand that the phase-field be approximately constant 
in the solid in spite of the periodic structure. If the solid structure is aligned with the 
computational domain in such a way that the global spatial averages are taken parallel to 
atomic layers, then the parameter e controlling the width of the average in the orthogonal 
direction must be large enough to smooth the gaps between the atomic layers. In the 
numerical simulations the orientations of the FCC lattice with respect to the solid-liquid 
interface, and hence the planes of averaging, are precisely such that averages are computed 
parallel to atomic planes, as illustrated in Figure [22] In the present case the distance to the 
nearest neighbours in the FCC-lattice is around 1.02; with rj on the form ( 16 ) the parameter 
e must be taken greater than 0.43 to ensure that r] decreases with at most a factor 1/2 in 
half the distance to the nearest neighbour, which seems a reasonable demand. Figure [23] 
presenting computed phase-fields based on local averages of the density and the potential 
energy using e = 0.45, shows that the smoothing parameter has to be larger than this to 
avoid oscillations in the solid part. The phase-fields based on e = 0.70 in Figure [24] do not 
show these oscillations on the length scale smaller than the distance between atom layers. 



For the method to be reasonable, the lower bound on e must not hide an interface width 
in the phase-field that is sharp even on the atomic scale. In addition to the computations 
with e = 1.0, the phase field has been computed for e = 0.45, 0.70, and 2.0. The computed 
phase-fields in the regions around the interfaces, for both orientation 1 and 2, are shown 
in Figure |25| The comparison shows that the interface width varies with the smoothing 
parameter. It would not, however, become infinitely sharp in the limit when e goes to zero, 
even if the lower bound on e were disregarded. This is clear from the results presented in 
Figure [26] where, in addition to the values of e above, a phase-field obtained with e = 0.05, 
violating the lower bound, is shown around one of the interfaces in 01. This value of the 
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(a) B 



(b) B 




(c) Some bj-.s 

Figure 20: The computed average diffusion matrix B, for e — 1.0, using the same configu- 
rations from simulation 02 as in Figure |8(d)| and Figure |17(d)[ is shown in (a) . The square 
root B of B, as defined in (20) is shown in (b). The individual columns in B are the dif- 
fusion coefficient functions, bj, in the stochastic differential equation for the coarse-grained 
phase- field m*. Some of these column vectors have been plotted as functions of the space 
variable Xi in (c). The support of each bj is centred around the grid point xj. 
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(a) FCC, e=2.0 (b) FCC, e=1.0 (c) FCC, e=0.55 
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(d) Liquid, e=2.0 (e) Liquid, e=1.0 (f) Liquid, e=0.55 



Figure 21: The average diffusion coefficient functions, 6(Aa;i) — mean + Aa;i)| have 

been computed for different values of e, with the mean taken over points ar[ in the interior 
of the soHd and the hquid domains, respectively. The configurations used are the same as 
in Figure [20] The difference between the solid and liquid parts is small. 
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Distance between atomic planes: D=(1/2) r 



Distance between atomic planes: D=(2/3) r 











1; . 1 

• 




1 • 




f • 1 


r J 




















(a) Orientation 1 
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Figure 22: The distance between two adjacent atom layers in a perfect FCC lattice 
is ^Jl/2ro in orientation 1 and ^2 /3 ro in orientation 2, where is the nearest neigh- 
bour distance. 
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Fi gure 23) Computed density, pioct 

and potential energy phase fields for simulations 01 

and 02 using e = 0.45. 
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smoothing parameter, and the corresponding moUifier cutoff, Rc = 6-0.05 = 0.3, is so small 
that the contribution to the phase-field of an individual atom in the FCC lattice is restricted 
to an interval extending less than halfway to the next atom layer in either direction. Still the 
change in the phase-field, from strong oscillations in the solid to decaying oscillations around 
the average in the pure liquid, occurs gradually on a length scale corresponding to at least 
several atom layers and thus several times the artificial smoothing introduced by e. Figure [26| 
also shows that the interface region of the phase-field obtained with e = 0.45, 0.70, and 1.0 is 
wider than the transition region of a step function, representing an infinitely sharp interface, 
smoothed by a convolution with the moUifier using the corresponding e. For e = 2.0 the 
interface is very close to that of a mollified step function in both width and profile. The 
interface width of the smoothed step function is proportional to e and it is expected that 
the same will hold for the phase-field, m, if the smoothing parameter is increased beyond 
the present range. 
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Figure 24: Computed density, pioa and potential energy phase fields for simulations 01 
and 02 using e = 0.70. 
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Figure 25: The moUifier, 77, in the definition of the phase field, m, depends on the model 
parameter e. The width of the averaging is proportional to e, as illustrated in (a) which 
shows ?7 for e = 0.45, 0.7, 1.0, 2.0. 

The phase field, m, in the interface regions has been computed from 174 configurations with 
the four e-values listed above. In (b) and (c) the configurations are taken from simulation 
01, and in (d) and (e) from simulation 02. In each case the time interval between two 
successive configurations is 2.5 • 10~^, corresponding to 5 • 10^ time steps. Though the 
interface width in the computed phase-fields varies with e, it is not proportional to e in this 
range. 
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Figure 26: For the phase field based on local contributions to the potential energy the 
transition from solid to liquid occurs on a length scale of at least several nearest neighbour 
distances for any choice of the smoothing parameter e. 

The four subfigures are based on the same configurations from simulation 01 as were 
Figure |25(b)f|25(c)l The oscillating curve present in all subfigures is the computed phase- 
field, m, using e = 0.05 with a cutoff of rj at 0.3. The nearest neighbour distance is 
approximately 1 and, for the present orientation of the FCC structure with respect to the 
a;i-axis, the xi-distance between the atomic layers becomes approximately 1 / \/2. Since 
the cutoff is less than half the distance between the atomic layers the phase-field would 
be exactly zero at the middle distance if the crystal were perfect and it is very close to 
zero here. The transition from the stable oscillation pattern in the solid to diminishing 
oscillations around the mean in the liquid is extended over a distance corresponding to at 
least four or five atomic layers in the solid. 

The phase-field, m, for e = 0.45, 0.70, 1.0, and 2.0 is shown as the heavy solid curve in 
subfigures (a)-(d). For reference the convolutions f{y)ri{x — y) dy of a sharp interface, 
given by the step function f{y) = miiqlR- (y) — r7iFcclR+(y), and the moUifier using the 
respective e-value is included as the heavy dashed curve. For the smaller e-values the 
mollified step function is significantly sharper than the corresponding phase-field. 
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(c) Rescaling factors — the accuracy is approximately ±0.05. 



Figure 27: For an interface given by the convolution of a sharp step function and the 
moUifier, as in Figure |26| the interface width is directly proportional to e, since interface 
profiles corresponding to different e are identical up to affine coordinate transformations 
around the interface, Xif, that is: (f)e2i^{x — Xif) + xn) = 4>ei{x). On a sufficiently large 
scale the same scaling of the interface width can be expected from the phase-field m obtained 
from MD simulations. This is not the case when e is of the order of the nearest neighbour 
distance; then the interface width grows more slowly than e. One way to quantify this 
statement is to consider the tabulated phase-field, mj^, using the parameter value ei, as 
given data to be approximated by the phase-field, m^^, based on the parameter value £2! 
the allowed approximations use affine coordinate transformations y{x) — Ci{x — Cq) + cq oi 
the independent coordinate. The data points ((x^), m^j^ (x^)) are taken from the interior 
of an interface, msoiid < m-o < '^^^{xk) < mi < muquid, and the function m^^ is defined 
by linear interpolation between tabulated values. A least squares approximation of the 
overdetermined system ■m^r^{y{xk)) = me^{xk) for cq and ci gives a value of the scaling 
factor ci to be compared to £2/61. 

Subfigures (a) and (b) show two examples for the interface in Figure |25(b)| The circles, 
o, denote the reference data points, the solid line shows the linear interpolation of the 
tabulated values for the approximating phase-field, and the line marked with crosses, x , is 
the least square approximation. 

The table (c) shows the scaling constants obtained after averaging over all four interfaces 



in Figure 25(b) ■25(e) The corresponding quotients ^ are included in parenthesis for 
reference. 
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A Explicit Calculation of Drift and Diffusion Functions 



Let the total potential energy be 

JV 
i=l 

where 

N 



For the phase-field 

m(x; X) = ^mi{X)r]{x - XJ, 



2 



where the particle positions X e M.^^ solve the Ito stochastic differential equation 

dX* = -VxU{X*) dt+ y/2kBT dW\ 

Ito's formula gives 

JV Af 3 

j=i j=i fe=i 

with 

aj{x; X) = -Vx,TO(a;; X) ■ Vx,U{X) + fcBTVx, • Vx,m(a;; X) (26) 

and 

= y2fcBTVx,m(x;X). (27) 

Introducing the total force, Fj, acting on particle j, and the contributions from individual 
pairs, fij, 

N 

F^{X) = -Vx,U{X)= J2 M^)^ 

/,,(X) = $'(||X,-X,||)^1^^, 
the gradient of rrii with respect to the position of particle j is 

AT 



Vx,miiX) = ^ ^ \/xM\\^^-^k 

1 ^ 1 

= <5..^ E Vx,*(||X,--X,||) + (l-(5,,)-Vx,1>(||X,-X^. 



= -5,,-^F,(X)-(l-(5,,)^/,,(X), 
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where 5ij is the Kronecker delta: 6ij = i = j,6ij = 0,if i 7^ j. The gradient of the 
phase-field variable with respect to the position of particle j is 



TV 



= -mj{X)V^rj{x- X^) 

N N 

- ^Y.^^JF,{XUx - X^) -^(1 - S,,)f,,{X)7^{x - X^) 

1=1 1=1 
= -V^{mj{X)'n{x-Xj)) 

1 1 N 

- F,{XHx^X^)-- ^ f,,{X)r,{x~X^). 



Introducing the notation —Gj for the divergence of the force Fj with respect to Xj and the 
notation gij for the individual contributions, 

N N 

G,{X) = -Vx, ■ F,{X) = - ^ Vx, ■ MX) = J2 

g,,{X) = <^"{\\X, - X^W) + ^'{\\X, - X^-ll)^^-^, 

the divergence of gradient of phase field variable with respect to the position of particle j 
becomes 



Vx, •Vx,m(x; X) = -Vx, • {m^{X)V^ri{x - X^)) 

1 1 ^ 

- -Vx, • (,F,{X)ij{x - X^)) -- J2 {.UjiXMx - X,)) 



Vx.m^iX) ■ V,7j{x ~ X^) - m,(X)Vx, • V,7y(a: - X^) 
Vx, • F,{X)rj{x - X^) - \f,{X) ■ Vx,77(x - X^) 



2 2 

1 ^ 

2 E Vx, •/.,(^)r7(a;-X,) 



1 
2 

+ ^G,(X)r7(x - X^.) + ^F,iX) ■ V,77(x - X^) 

N 



+ 1 J2 9v{XHx~X,) 

V,-V, {mj{X)r^{x - X^)) + V, • {Fj{X)t^{x - X^)) 

1 1 ^ 

+ -G,(X),7(x-X^) + - 5] g.,,{X)r,{x-X,). 



Using the explicit expressions for \7 Xjin-ix; X) and \7 Xj'^ Xj'nT'{x; X), the components (26l 
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of the drift become 

a,{x;X) = V,{m^{XUx - X^)) ■ (-F.iX)) + ^FjiX^x - X^) ■ 



= -V, • im,{X)F,{X)r,{x - X^)) - -\\Fj{X)\\'rj{x - X^) 



+ ^ E fijiXHx - X,) . i-F^iX)) 

+ kBTVx,-"^x,m{x;X) 

1 ^ 

-2 E MX) ■ FjiXMx - X,) 

+ kBTV^-V^ {mj{X)T]{x - Xj)) + knTV^ ■ {Fj{X)r,{x - Xj)) 

1 1 ^ 

+ kBT-Gj{X)rj{x-X^)+kBT- ^ ff,,(X)r?(a; - X,). 

^ kBTV,-V^{mj{X)7j{x - X^)) 
+ V, • ({kBT - mj{X))Fj{X)r]{x - Xj) 

+ ^-{kBTG,iX)-\\FjiX)\\').,y.-..^, 

N 



-{kBTG,{X)-\\FjiX)\\^)rjix-X^) 
1 ^ 

+ x E (kBTgijiX) - f,j{X) ■ Fj{X))r,{x - X,) 



2 , 



so that, after summing over j, 

a{x; X) = A;bTVx- V^m(a;; X)+V^- ai{x; X) + ao{x; X) 



with 



and 



N 

ai{x; X) = E(^bT - mj{X))Fj{X)r]{x - Xj) 



ao{x;X) = (^knTGjiX) - r?(x - X.) 

N N 

-2E E MX)-Fj{XUx-X^). 



Using the one-dimensional moUifier 



ry(a;) = r]{xi) = constant • exp ( — ^ ^— ) ) , (28) 
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that only varies in the xi -direction, the expression for the drift reduces to 

52 d 
a{x; X) — k^T—^m{x; X) + - — ai{x; X) + ao(x; X) 

with 

N 

a^{x- X) = J2ikBT - m,{X))[F,{X)lrj{x - X^), 
i=i 

where is the xi component of Fj{X). 



For the purpose of computing an approximation of 

rT N 3 



1 



it is not practical to postpone the differentiation of the moUifier with respect to the space 
varible. Using the choice (28 1, the gradient of the moUifier can be expressed in terms of the 
moUifier itself as 

V^x-Xj) = -^jj{x-Xj) (^[x~X^]^,0,Oy . 
Then the expression for W Xjrn{x; X) becomes 

Vx,m(x;X) = (^^^^ [[x - X^]^,0,oY~ ^^-'(^)) " 
1 ^ 

"2 E MXHx-X,) 



and, using the diffusion component (27 1, 

3 



k=l 



where 

P]{x,y;X) 



mj{X) 
2e2 

2e2 



and 



N 

[x^X^],rj{x-X^)([F,{X)],rj{y-X^)+ ^ [f,,{X)],Tj{y - X^) 

N 

[y ~ X^],7^{y - X^)([F,{X)lrj{x - X^) + ^ [f,,{X)lrj{x - X,) 



1 / N 

{x,y;X) = -(F,{X)r,{x^X^)+ ^ f,,{XMx - X,) 



N 



F,{X)rj{y~X^)+ ^ f.,,{XUy - X,) 
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